#include "headers.hh"
#include "SysUtils.h"

void cutoffID51(const TList * ll)
{
  const Int_t ngr = ll->GetEntries();
  for(Int_t igr=0; igr<ngr; igr++){
    TGraphAsymmErrors * tmpgr = (TGraphAsymmErrors *) ll->At(igr);
    Int_t istart = -999;
    for(Int_t ipt=0; ipt<tmpgr->GetN(); ipt++){
      //cut off at 0.6 gev
      if(tmpgr->GetX()[ipt]>0.6){
        istart=ipt;
        break;
      }
    }
    if(istart<0){
      printf("cutoffID51 bad istart!\n"); exit(1);
    }

    const Int_t newn = tmpgr->GetN()-istart;
    printf("cutoffID51 oldn %d newn %d\n", tmpgr->GetN(), newn);
    Int_t ipt=0;
    for(ipt=0; ipt<newn; ipt++){
      const Int_t oldi = ipt+istart;

      if(tmpgr->GetX()[oldi]>12)//high end electron statistics
        break;

      tmpgr->SetPoint(ipt, tmpgr->GetX()[oldi], tmpgr->GetY()[oldi]);
      tmpgr->SetPointError(ipt, tmpgr->GetEXlow()[oldi], tmpgr->GetEXhigh()[oldi], tmpgr->GetEYlow()[oldi], tmpgr->GetEYhigh()[oldi]);
      printf("cutoffID51 %s %d/%d %e %e\n", tmpgr->GetName(), ipt, newn, tmpgr->GetX()[oldi], tmpgr->GetY()[oldi]);
    }
    tmpgr->Set(ipt);
  }
}

TH2D * get2dHist(const TList * ll, const TString tag, const Bool_t *badgrs)
{
  TH2D * hh = new TH2D(Form("%shh",tag.Data()),"", 50, 0.1, 100,200,-1,1);
  SysUtils::BinLogX(hh->GetXaxis());
  hh->GetXaxis()->SetRangeUser(0.2, 40);

  const Int_t ngr = ll->GetEntries();
  for(Int_t igr=0; igr<ngr; igr++){
    if(badgrs && badgrs[igr])
      continue;

    TGraphAsymmErrors * tmpgr = (TGraphAsymmErrors *) ll->At(igr);
    for(Int_t ipt=0; ipt<tmpgr->GetN(); ipt++){
      hh->Fill(tmpgr->GetX()[ipt], tmpgr->GetY()[ipt]);
    }
  }

  return hh;
}

Int_t setCpbad(Bool_t badgrs[], const Int_t nfile, const Double_t cparray[], const Double_t nsigma)
{
  Double_t sum = 0;
  Double_t w2 = 0;
  Int_t ndata = 0;

  for(Int_t ii=0; ii<nfile; ii++){
    if(badgrs[ii]){
      continue;
    }

    sum += cparray[ii];
    w2 += cparray[ii]*cparray[ii];
    ndata++;
  }
  
  if(ndata<2){
    printf("setCpbad ndata <2!! %d\n", ndata); 
    //test
    //exit(1);
    return -1;
  }

  const Double_t mm = sum/ndata;
  //if rms>5 per mil the distribution is very bad!
  const Double_t rms = TMath::Min(0.005, TMath::Sqrt(w2/ndata-mm*mm));
  Int_t nbad=0;
  for(Int_t ii=0; ii<nfile; ii++){
    if(badgrs[ii]){
      continue;
    }

    if(TMath::Abs(cparray[ii]-mm)>nsigma*rms){
      badgrs[ii]=kTRUE;
      nbad++;
    }
  }

  return nbad;
}

Int_t setbadid(const TList * ll, Bool_t badgrs[], TList *lout,  const Int_t dpt, const Double_t nsigma)
{
  //npt should be the smallest in the list because the xbinning also changes and npt is not constant!!
  Int_t npt=10000;
  for(Int_t igr=0; igr<ll->GetEntries(); igr++){
    const Int_t tt = ((TGraphAsymmErrors *) ll->At(igr))->GetN();
    if(tt<npt){
      npt=tt;
    }
  }

  //=================

  Int_t nbad = 0;

  static Int_t counter=0;

  //do it quasi bin-by-bin
  for(Int_t ipt = 0; ipt<npt; ipt+=dpt){

  //=======================

  const Int_t ngr = ll->GetEntries();

  Double_t pars[ngr];
  Int_t npar[ngr];


  TH1D * hgr = new TH1D(Form("hgr%dipt%d", counter, ipt),"",(Int_t) pow(2.,11.), 0,1); lout->Add(hgr);
  counter++;

  for(Int_t igr=0; igr<ngr; igr++){
    pars[igr]=0;
    npar[igr]=0;

    if(badgrs[igr])//only consider good ones
      continue;

    TGraphAsymmErrors * tmpgr = (TGraphAsymmErrors *) ll->At(igr);
    for(Int_t ik=0; ik<dpt; ik++){
      if(ipt+ik>=npt)
        break;

      pars[igr]+=tmpgr->GetY()[ipt+ik];
      npar[igr]++;
      //printf("test igr %d ipt %d ik %d -- %d %d %e %d\n", igr, ipt, ik, ipt+ik, npt, tmpgr->GetY()[ipt+ik],  npar[igr]);

    }

    if(npar[igr]==0){
      printf("setbadid npar=0! igr %d npt %d counter %d\n", igr, npt, counter); 
      exit(1);
    }
    pars[igr]/=npar[igr];
    hgr->Fill(pars[igr]);
    printf("setbadid %s counter %d igr %d pars %e\n", tmpgr->GetName(),   counter, igr, pars[igr]);
  }

  const Double_t mean = hgr->GetMean();
  const Double_t rms = TMath::Min(0.1, hgr->GetRMS());
  printf("setbadid mean %e rms %e nsigma %e\n", mean, rms, nsigma);
  for(Int_t igr=0; igr<ngr; igr++){
    if(!badgrs[igr]){
      if(TMath::Abs(pars[igr]-mean)>nsigma*rms){
        printf("setbadid counter %d ipt %d dpt %d igr %d setbad pars %e mean %e nsig %e!\n", counter, ipt, dpt, igr, pars[igr], mean, TMath::Abs(pars[igr]-mean)/rms);
        badgrs[igr]=kTRUE;
        nbad++;
      }
    }
  }

  }

  return nbad;

  /*
  const Double_t pthres = 10;
  const Double_t ythres = 0.1;
  for(Int_t igr=0; igr<ngr; igr++){
    TGraphAsymmErrors * tmpgr = (TGraphAsymmErrors *) ll->At(igr);
    for(Int_t ipt=0; ipt<tmpgr->GetN(); ipt++){
      if(tmpgr->GetX()[ipt]>pthres && tmpgr->GetY()[ipt]<ythres){
        badgrs[igr]=kTRUE;
      }
    }
  }
  */

  /*
  TH1D * hrec=0x0;
  Int_t clsid = -999;
  Double_t cut = -999;
  TH1D * hout = SysUtils::GetMajorCluster(hkmean, hrec, clsid, cut);
  lout->Add(hrec);
  lout->Add(hout);

  //for kaon, only choose clsid==2
  if(clsid==2){
    for(Int_t ii=0; ii<ngr; ii++){
      if(meankaon[ii]<cut){
        badgrs[ii]=kTRUE;
      }
    }
  }
  */

  /*
  Int_t ids[ngr];
  TMath::Sort(ngr, meankaon,  ids, kFALSE);

  const Double_t frac = 0;//0.20; //~5/129; //0.5; too big, sys. error will be too small
  const Int_t nbad = ngr*frac;
  for(Int_t ii=0; ii<nbad; ii++){
    badgrs[ids[ii]]=kTRUE;
    //printf("test0 %d %d bad\n", ii, ids[ii]);
  }
  */
}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=5){
    printf("argv!=5 \n"); return 1;
  }

  ifstream input(argv[1]);
  if(!input){
    printf("getGrSys file not found \n%s\n\n",argv[1]); exit(1);
  }

  const Int_t grid = atoi(argv[2]);
  const Int_t ksetbad = atoi(argv[3]);
  const Int_t kCheckCp  = atoi(argv[4]);

  const Int_t ntype = 4;

  const Int_t nfmax = 2000;
  TString fins[nfmax];
  TString buff;
  Int_t nfile=0;
  while(input.good()){
    input>>buff;
    if(buff!="" && !buff.Contains("#")){
      fins[nfile]=buff;

      printf("reading %d %s\n", nfile, fins[nfile].Data());

      nfile++; 

      if(nfile>nfmax){
        printf("getGrSys nfile>nfmax %d %d\n", nfile, nfmax); return 1;
      }   
    }
  }

  TList *ll[ntype];
  TList *lid10[ntype];
  for(Int_t itype=0; itype<ntype; itype++){
    ll[itype]=new TList;
    lid10[itype]=new TList;
  }

  Double_t cparray[nfile];

  for(Int_t ifile=0; ifile<nfile; ifile++){
    printf("ifile: %d %s\n", ifile, fins[ifile].Data());

    cparray[ifile]=-999;

    TFile *ff=new TFile(fins[ifile]);
    if(!ff){
      printf("getGrSys no ff!! %d %s\n", ifile, fins[ifile].Data()); exit(1);
    }

    for(Int_t itype=0; itype<ntype; itype++){
      printf("itype: %d\n", itype);

      const TString gname(Form("grmn2dp%d%02d",itype, grid));
      TGraphAsymmErrors *gr = (TGraphAsymmErrors*) ff->Get(gname);
      if(!gr){
        printf("getGrSys gr not found!! %d %d %s %s\n", itype, ifile, fins[ifile].Data(), gname.Data());exit(1);
      }
      gr->SetName(gname+Form("f%d",ifile));
      ll[itype]->Add(gr);

      //00: is p base fraction
      //10: can be rebinned and therefore not suitable to make removal
      //-> all wrong!
      //p base has different binning and therefore no way to cut on the same p!!
      //must use pT
      const TString tmpgn(Form("grmn2dp%d%02d", itype, 10));
      TGraphAsymmErrors *tmpg10= (TGraphAsymmErrors*) ff->Get(tmpgn);
      if(!tmpg10){
        printf("no tmpg10 for setbadid!! should not matter if no setbadid\n");
      }
      else{
        tmpg10->SetName(tmpgn+Form("f%d", ifile));
        lid10[itype]->Add(tmpg10);
      }
    }

    //------------
    if(kCheckCp){
      TString scpf = fins[ifile];
      scpf.ReplaceAll("execlocal/parmn2d_","outmn2d_post_");
      scpf.ReplaceAll("execqsub/parmn2d_","outmn2d_post_");
      scpf.ReplaceAll("execlocal/badfile1parmn2d_","outmn2d_post_");
      scpf.ReplaceAll("execqsub/badfile1parmn2d_","outmn2d_post_");
      scpf.ReplaceAll("execlocal/badfile2parmn2d_","outmn2d_post_");
      scpf.ReplaceAll("execqsub/badfile2parmn2d_","outmn2d_post_");
      scpf.ReplaceAll("root","#root");
      scpf=scpf(0, scpf.First("#"));
      scpf+="_save.root";
      TFile * fcp = new TFile(scpf);
      if(!fcp->IsOpen()){
        printf("no fcp %d %s %s\n", ifile, fins[ifile].Data(), fcp->GetName());exit(1);
      }
      TTree *tcp = (TTree*) fcp->Get("tree");
      if(!tcp){
        printf("no tree in fcp %s\n", fcp->GetName()); fcp->ls(); exit(1);
      }
      tcp->SetBranchAddress("bb7", &(cparray[ifile]));
      tcp->GetEntry(0);
      printf("C_p %d: %e\n", ifile, cparray[ifile]);
    }
  }

  printf("getGrSys nfile %d\n", nfile);
  TList *lout=new TList;

  //===============================================================
  if(grid==51){
    for(Int_t itype=0; itype<4; itype++){
      cutoffID51(ll[itype]);
    }
  }
  //===============================================================

  const Int_t ngr = nfile;
  Bool_t badids[ngr];
  for(Int_t ii=0; ii<ngr; ii++){
    badids[ii]=kFALSE;
  }

  //all grid should have the same bad list
  //pi, K, p, e
  //pi only!
  if(ksetbad && nfile>1){
  const Int_t neworder[]={1,2,0,3};//,2,0,3};
  const Int_t nt=sizeof(neworder)/sizeof(Int_t);
  for(Int_t io=0; io<nt; io++){
    const Int_t itype = neworder[io];

    Int_t nfinalbad=0;
    const Int_t niter = 50;
    for(Int_t ii=0; ii<niter; ii++){
      nfinalbad=0;

      //bin-by-bin for dpt = 1, 3-sigma
      printf("badtest io %d ii %d %d\n",io, ii,  itype);
      const Int_t dp = 5;
      const Double_t nsigma = 3;
      nfinalbad+= setbadid(lid10[itype], badids, lout, dp, nsigma);

      //Int_t setCpbad(Bool_t badgrs[], const Int_t nfile, const Double_t cparray[], const Double_t nsigma)
      if(kCheckCp){
        const Double_t outCpbad = setCpbad(badids, nfile, cparray, nsigma);
        if(outCpbad<0){
          TFile * fout=new TFile("check.root","recreate");
          lout->Write();
          fout->Save();
          fout->Close();
          exit(1);
        }
        nfinalbad+= outCpbad;
      }

      if(nfinalbad==0){
        break;
      }
    }
    if(nfinalbad){
      printf("getGrSys still find bad in final iteration %d! %d\n", nfinalbad, niter); exit(1);
    }
  }
  }

  Int_t nacc = 0;
  TH1D * hcpbefore=new TH1D("hcpbefore","",200,-0.06,0.06); lout->Add(hcpbefore);
  TH1D * hcpafter=new TH1D("hcpafter","",200,-0.06,0.06); lout->Add(hcpafter);
  for(Int_t ii=0; ii<ngr; ii++){
    if(!badids[ii]){
      printf("getGrSys accepted file %d: id %d %s\n", nacc, ii, fins[ii].Data());
      hcpafter->Fill(cparray[ii]);
      nacc++;
    }
    hcpbefore->Fill(cparray[ii]);
  }
  printf("getGrSys %s acceptedfraction %e nacc %d nfile %d\n", gSystem->pwd(), (Double_t)nacc/(Double_t)nfile, nacc, nfile);

  if(grid==10 || grid==60 || grid==40 || grid==90 || grid == 1 || grid ==2){
    for(Int_t itype=0; itype<ntype; itype++){
      TH2D * hdummy = 0x0;
      hdummy = get2dHist(ll[itype],    Form("ptbefore%d",itype),  0x0); lout->Add(hdummy);
      hdummy = get2dHist(ll[itype],    Form("ptafter%d",itype),  badids); lout->Add(hdummy);
    }
  }
  //===============================================================

  //===================================================================
  Bool_t kfit = kFALSE;
  for(Int_t itype=0; itype<ntype; itype++){
    /*
    if(grid/10<5 && itype<3){
      kfit = kTRUE;
    }
    else {
      kfit = kFALSE;
    }
    */

    printf("getGrSys kfit %d grid %d itype %d\n", kfit, grid, itype);
    const Double_t hmin = 0; 
    Double_t hmax = 1;
    if(grid%10==1){//01
      hmax = 200;
    }
    else if(grid%10==2){//02
      hmax = 0.1;
    }

    SysUtils::GetGrSysMode0(ll[itype], Form("%d",itype), badids, kfit, hmin, hmax);
  }

  //======================================================================

  TGraphAsymmErrors * grsys[ntype];
  TGraphAsymmErrors * grstat[ntype];
  for(Int_t itype=0; itype<ntype; itype++){
    grsys[itype] = (TGraphAsymmErrors*) ll[itype]->FindObject(Form("grsys%d", itype));
    grstat[itype] = (TGraphAsymmErrors*) ll[itype]->FindObject(Form("grstat%d", itype));
    if(!grsys[itype] || !grstat[itype]){
      printf("no grsys[itype] %d\n", itype); return 1;
      ll[itype]->ls();
    }
  }

  //check consistency of npoint between different itype
  Int_t npt = -999;
  for(Int_t itype=0; itype<ntype; itype++){
    if(grid%10==8 && itype==3)
      continue;

    if(grid%10==1 || grid%10==2)
      continue;

    if(npt<0){
      npt = grsys[itype]->GetN();
    }
    else if(npt!=grsys[itype]->GetN()){
      printf("getGrsys bad npt! itype %d npt %d getn %d\n", itype, npt, grsys[itype]->GetN()); return 1;
    }
  }

  for(Int_t itype=0; itype<ntype; itype++){ 
    for(Int_t ip=0; ip<npt; ip++){      
      printf("getGrSys finalgraph itype %d ip %d x %e exl %e exh %e y %e stat eyl %e eyh %e sys eyl %e eyh %e\n", itype, ip, grsys[itype]->GetX()[ip], grsys[itype]->GetEXlow()[ip], grsys[itype]->GetEXhigh()[ip], grstat[itype]->GetY()[ip], grstat[itype]->GetEYlow()[ip], grstat[itype]->GetEYhigh()[ip], grsys[itype]->GetEYlow()[ip], grsys[itype]->GetEYhigh()[ip]);
    }
  }

  //check sum = 1
  //not needed if kfit any way is not used
  /*
  if(grid%10==0){
    for(Int_t ip=0; ip<npt; ip++){
      Double_t suma = 0;
      Double_t tmpfrac[]={-999,-999,-999,-999};
      for(Int_t itype=0; itype<ntype; itype++){
        tmpfrac[itype] = grsys[itype]->GetY()[ip];
        suma+= tmpfrac[itype];
      }
      if(fabs(suma-1)>1e-3){
        printf("getGrSys error suma!=1 ip %d %e , %e %e %e %e\n", ip, suma-1, tmpfrac[0], tmpfrac[1], tmpfrac[2], tmpfrac[3]); //test exit(1);
      }
    }
  }
  */
  //======================================================================

  TFile *fout=new TFile(Form("outgrsys_id%02d.root", grid),"recreate");
  for(Int_t itype=0; itype<ntype; itype++){
    //ll[itype]->Write(Form("ll%d",itype), TList::kSingleKey);
    ll[itype]->Write();
  }
  lout->Write();
  fout->Save();
  fout->Close();

  printf("getGrSysdone!\n");

  return 0;
}

